The real power of GrADS lies in its data analysis capabilities. These capabilites are accessed through GrADS internal or user defined external functions (so called "udfs") applied to an expression (GrADS talk for a grid of data, e.g, var.1(t=1)).
These functions are biult into GrADS.
Math: abs, acos, asin, atan2, cos, exp, log, log10, pow, sin, sqrt, tan
Filtering: smth9
Station Data: gr2stn, oacres, stnave, stnmin, stnmax
Finite Difference: cdiff
Special Purpose: tloop
Meteorological: tvrh2q, tvrh2t
Functions are invoked by name, with their arguments separated by commas and enclosed in parentheses.
Expressions are typically one of the arguments supplied to a function. Functions may be nested. Some functions modify the dimension environment when they operate.
Users may write their own GrADS functions in the computer language of their choice, and have them available from the GrADS expression facility (via the display command). This feature is currently in testing.
aave(expr,xdim1,xdim2,ydim1,ydim2)
Takes an area average over an X-Y region. Usually better than using nested ave functions.
expr any valid GrADS expression.
xdim1 starting dimension expression for the X dimension.
xdim2 ending dimension expression for the X dimension.
ydim1 starting dimension expression for the Y dimension.
ydim2 ending dimension expression for the Y dimension.
Usage Notes
1) In the absence of missing data values, aave gives the same result as nested ave functions in the X and Y dimensions:
ave(ave(expr,x=1,x=72),y=1,y=46)
being equivalent to:
aave(expr,x=1,x=72,y=1,y=46)
in terms of the numerical result. The aave function is more efficient.
2) When there are missing data values, the aave function does not return the same result as nested ave functions. To see this, consider the small grid:
6 18 3 5
10 10 10 10
12 U U U
where U represents the missing data value. If we apply nested ave functions, the inner ave will provide row averages of 8, 10, and 12. When the outside ave is applied, the result will be an average of 10. When aave is used, all the values participate equally (in this case, we are assuming no weights applied to the final average), and the result is 84/9 or about 9.33.
3) The aave function assumes that the world coordinates are longitude in the X dimension and latitude in the Y dimension, and does weighting in the latitude dimension by the delta of the sin of the latitudes. Weighting is also performed appropriately for unequally spaced grids.
4) The aave function always does its average to the exact boundaries specified, in world coordinates. This is somewhat different from the ave function, where the -b flag is used to get this behavior. If the boundaries specified via the dimension expressions do not fall on grid boundaries, then the boundary values are weighted appropriately in the average.
When grid coordinates are used in the dimensions expressions, then they are converted to
world coordinates for the boundary to be determined. This conversion is done using the scaling of the default file. Note that the conversion is done using the outside grid boundary, rather than the grid center. For example:
aave(expr,x=1,x=72,y=1,y=46)
Here the boundary would be determined by using the grid values 0.5, 72.5, 0.5, and 46.5 which would be converted to world coordinates. If we assume that x=1 is 0 degrees longitude and x=72 is 355 degrees longitude, then the averaging boundary would be -2.5 to 357.5 degrees, which would cover the earth. In the Y dimension, when the boundary is beyond the pole, the aave function recognizes this and weights appropriately.
Examples
1) See the tloop function for an example of creating a time series of area averages.
2) An example of taking an area average of data only over land, given a mask grid:
aave(maskout(p,mask.3(t=1)),x=1,x=72,y=1,y=46)
In this case, it is assumed the mask grid has negative values at ocean points.
abs(expr)
Takes the absolute value of the result of expr. Operates on both gridded and station data. Missing data values do not participate.
acos(expr)
Applies the cos-1 function to the result of expr. Values from expr that exceed 1 or are less than -1 are set to missing. The result of the acos function is in radians.
asin(expr)
Applies the sin-1 function to the result of expr. Values of expr that exceed 1 or are less than -1 are
set to missing in the final result. The result of the asin function is in radians.
atan2(expr1,expr2)
Applies the tan-1 function to the result of the two expressions, using the formula:
where y is expr1 and x is expr2. Situations where x and y are both zero are valid; the result is arbitrarily set to zero. The result of the atan function is in radians.
ave(expr,dexpr1,dexpr2<,tincr<,flags>>)
Averages the result of expr over the specified range of dimensions starting at dexpr1 and ending at dexpr2. If the averaging dimension is time, an optional time increment tincr may be specified.
expr is any valid GrADS expression.
dexpr1 is the start point for the average, specified as a standard GrADS dimension expression.
dexpr2 is the end point for the average. The dimensions of dexpr1 and dexpr2 must match.
tincr is a time increment for the average, if dexpr1 and dexpr2 are time dimension.
flags The following flags are valid:
-b The boundry flag indicates the average should be taken to the exact boundaries specified in dexpr1 and dexpr2, rather than nearest grid points. See Usage Notes below.
Usage Notes
1) The limits and interval over which to take the average are determined by the scaling of the
default file. Conversions of dexpr1 and dexpr2 to grid coordinates are performed using
the scaling of the default file. See the Examples below for an example of what this means.
2) The average is weighted for non-linear grid intervals. Averages over latitude are weighted
by the delta of the sine of the latitudes at the edge of the grid box. The edges of the grid box are always defined as being the midpoint between adjacent grid points.
3) If dexpr1 and dexpr2 are specified in world coordinates, the coordinates are converted to
the nearest integer grid coordinates based on the scaling of thE default file. The average is then performed over the range of these grid coordinates. The end points are given normal weighting, unless the -b flag is specified.
Examples
All the examples use the two example descriptor files given at the beginning of this section. For the following examples, the dimension environment is X-Y varying; Z-T are fixed.
1) Consider the following average, when the default file is file #1:
ave(z.2,t=1,t=10)
We are averaging a variable from file #2, but using the scaling from file #1. File #1 has a time interval of 6 hours, but file #2 has a time interval of 12 hours. The average will thus attempt to access data from file #2 for times that are not available, and an error will ocurr. To avoid this, the default file should be set to file #2: set dfile 2
2) The average:
ave(z,t=1,t=120,4)
will average only 00Z reports from file #1, since the time increment is 4, which for this file is 24 hours.
3) If you attempt to take a zonal average as follows:
ave(z,lon=0,lon=360)
the world coordinates will be converted to grid coordinates, here X varying from 1 to 181, and the grid point at longitude 0 (and 360) will be used twice in the average. To have the end points of this average weighted properly, use the -b flag:
ave(z,lon=0,lon=360,-b)
or average using the grid coordinates directly:
ave(z,x=1,x=180)
4) You can nest averaging operations:
ave(ave(z,x=1,x=180),y=1,y=46)
In this case, to take an areal average. Note that for areal averaging, the aave function is better. See the aave function description.
When nesting averages, the order of the nesting can have a dramatic affect on performance.
Keep in mind the ordering of the data in a GrADS file: X varies the fastest, then Y, then Z, then T. When nesting averages, put the faster varying dimension within the inner average:
set lon -90
set lat -90 90
set lev 1000 100
d ave(ave(t,x=1,x=180),t=1,t=20)
This average would be more efficient than, for example:
ave(ave(t,t=1,t=20),x=1,x=180)
although the final numerical result would be the same.
5) The use of the define command can make certain operations much more efficient. If you want to calculate standard deviation, for example:
sqrt(ave(pow(ave(z,t=1,t=20)-z,2),t=1,t=20))
would be correct, but the inside average would be calculated 20 times. Defining the inside average in advance will be substantially faster:
define zave = ave(z,t=1,t=20)
d sqrt(ave(pow(zave-z,2),t=1,t=20))
cdiff(expr,dim)
Performs a centered difference operation on expr in the direction specified by dim. The difference is done in the grid space, and no adjustment is performed for unequally spaced grids. The result value at each grid point is the value at the grid point plus one minus the value at the grid point minus one. The dim argument specifies the dimension over which the difference is to be taken, and is a single character: X, Y, Z, or T.
Result values at the grid boundaries are set to missing.
Examples
1) The cdiff function may be used to duplicate the calculation done by the hcurl function:
define dv = cdiff(v,x)
define dx = cdiff(lon,x)*3.1416/180
define du = cdiff(u*cos(lat*3.1416/180),y)
define dy = cdiff(lat,y)*3.1416/180
display (dv/dx-du/dy)/(6.37e6*cos(lat*3.1416/180))
The above example assumes an X-Y varying dimension environment. Note that the intrinsic variables lat and lon give results in degrees and must be converted to radians in the calaculation. Also note the radius of the earth is assumed to be 6.37e6 meters thus the U and V winds are assumed to have units of m/sec.
2) Temperature advection can be calculated using the cdiff function as follows:
define dtx = cdiff(t,x)
define dty = cdiff(t,y)
define dx = cdiff(lon,x)*3.1416/180
define dy = cdiff(lat,y)*3.1416/180
display -1*( (u*dtx)/(cos(lat*3.1416/180)*dx) + v*dty/dy )/6.37e6
where the variable t is temperature, u is the U component of the wind, and v is the V component of the wind.
const(expr,constant<,flag>)
Some or all of the values in expr are set to the constant value, depending on the value of flag:
expr A valid GrADS expression.
constant A constant, given as an integer or floating point value. The value will be treated as as floating point.
flag If no flag is specified, all the non-missing data values in expr are set to the constant value to form the result. Missing data values are preserved as missing.
Other flag values:
-u All missing data values are set to the constant value. Non-missing data remains unchanged.
-a All possible data values in the result are set to the constant value, including missing data values.
Usage Notes
1) This functions operates on both gridded and station data.
Examples
1) The const function may be used to assign a new value to missing data, so that missing data may participate in operations:
const(z,0.0,-u)
2) The const function is useful when displaying plots using the linefill graphics output type when one of the lines needs to be a straight horizontal line:
set lon -90
set lat -90 90
set gxout linefill
set lev 500
display const(t,-20);t-273.16
3) The const function may be used to calculate the fraction of the globe convered by some value of interest. In this case, the portion of the globe covered by precipitation greater than 10mm/day is calculated as a time series:
set lon 0 360
set lat -90 90
set t 1 last
define ones = const(const(maskout(p,p-10),1),0,-u)
set x 1
set y 1
display tloop(aave(ones,lon=0,lon=360,lat=0,lat=360))
Note that we first create a defined array that contains 1 wherever the precip value is greater
than 10, and 0 whever the precip value is less than 10. This is done via nested functions, where we first use the maskout function to set all values less than 10 to missing. We then use a const function with no arguments to set all non-missing values to 1, then use a const function with the -u flag to set all the missing data values to 0. The aave function is used to calculate an area weighted average. Since we are averaging zeros and ones, the result is the fraction of the area where there are ones. See the tloop function for a description of how to perform time series of areal averages.
cos(expr)
Takes the cosine of the expr. Values are assumed to be in radians. Works on both gridded and station data.
exp(expr)
Performs the e**x operation, where expr is x. Works on both gridded and station data.
gr2stn(grid_expr,stn_expr)
Performs an interpolation from grid space back to station locations, where:
grid_expr is any valid GrADS expression that gives a grid result. The interpolation will be done using this gridded data.
stn_expr is any valid GrADS expression that gives a station data result. The interpolation will be done to these station locations, the value of the station data is not used.
The result of the function is station data. The interpolation is done bi-linearly within the grid space. No weighting is done to account for real-world coordinate systems.
Examples
1) To examine the difference between an analysis (ie, gridded data) and the original observations, one could:
d t.3-gr2stn(t.1,t.3)
where file 1 is gridded data, and file 3 is station data. The result would displayas differences at the station locations.
2) If one wanted to display the difference calculated in Example 1 as a contour field, one can use the oacres function to do a quick analysis of the station values:
d oacres(t.1,t.3-gr2stn(t.1,t.3))
hcurl(uexpr,vexpr)
Calculates the vertical component of the curl (ie, vorticity) at each grid point using finite differencing on the grids provided. It is assumed that uexpr gives the U Wind component, and that vexpr provides the V Wind component.
Usage Notes
1) The alghorithm used for the finite difference calculation is described as an Example for the cdiff function.
2) The function assumes an X-Y varying dimension environment, and will not operate unless that is the case. The define command can be used in conjunction with the hcurl function to create 3 or 4 dimensional fields of vorticity, from which vertical cross-sections could be displayed.
3) The boundaries of the grid are set to missing.
4) The radius of the earth used in the calculation is in meters; thus the units of the wind expressions provided would normally be m/s.
Examples
1) To display the vorticity:
d hcurl(u,v)
2) If you want to display a vertical cross section of vorticity, you first need to calculate it over a 3-Dimensional region:
set lon 0 360
set lat -90 90
set lev 1000 100
define vort = hcurl(u,v)
set lon -90
display vort
hdivg(uexpr,vexpr)
Calculates the horizontal divergence using finite differencing. Exactly the same as hcurl in all other respects; see the Usage Notes and Examples above.
Usage Notes
1) The numerical stability of calculating horizontal divergence using finite differencing is very low. Please use the function with caution.
log(expr)
Takes the natural logarithm of the expression. May be used with gridded or station data. Values less than or equal to zero are set to missing in the result.
log10(expr)
Takes the logarithm base 10 of the expression. May be used with gridded or station data. Values less than or equal to zero are set to missing in the result.
mag(uexpr,vexpr)
Performs the calculation: sqrt(uexpr*uexpr+vexpr*vexpr). May be used with gridded or station data.
maskout(expr,mask)
Wherever the mask values are less than zero, the values in expr are set to the missing data value.
Works with gridded or station data. Where mask values are positive, the expr values are not modified. Thus the result of maskout is data with a possibly increased number of missing data values. The maskout function, in spite of its apparant simplicity, is extremely useful.
Examples
1) See the Examples for the const function for a description of using maskout to calculate the percentage of the globe covered by precipitation.
2) The maskout function can be used to cause part of the data to be ignored while doing another calculation. For example, if we have a land-sea mask, where sea values are negative, and we want to take some areal average of a quantity only over land:
d aave(maskout(p,mask.2),lon=0,lon=360,lat=0,lat=90)
3) People frequently have trouble using a mask grid, because it is often put into a seperate file, and given some arbitrary date/time and level. Thus, it is often necessary to locally override the dimension environment while using the mask grid:
d aave(maskout(p,mask.2(t=1)),lon=0,lon=360,lat=0,lat=90)
would probably be how Example 2 would have to be expressed in order to work, with the local override of t=1 specified on the mask data. See the documentation on how GrADS evaluates expressions within the dimension environment for more information.
oacres(grid_expr,stn_expr<,radii>)
A Cressman objective analysis is performed on the station data to arrive at a gridded result that represents the station data:
grid_expr An expression that has a gridded data result. The actual values of this result are ignored; the grid is used as a template to perform the analysis. The scaling of this grid must be linear in lat-lon.
stn_expr An expression that has a station data result. The station data is analyzed to the grid.
radii Optional radii of influence. Multiple radii are usually provided, seperated by commas. If not provided, default radii are used, in grid units: 10,7,4,2,1. The third radius specified is special, in that any grid points that do not have stations within that radius are set to the missing data value. See below for a discussion of the radii of influence.
The Cressman Analysis scheme is described in a paper in Monthly Weather Review, 1959. In summary, multiple passes are made through the grid at subsequently lower radii of influence. At each pass, a new value is determined for each grid point by arriving at a correction factor for that grid point. This correction factor is determined by looking at each station within the radius of influence from the grid point. For each such station, an error is determined as the difference between the station value and a value arrived by interpolation from the grid to that station. A distance weighted formula is then applied to all such errors within the radius of influence of the grid point to arrive at a correction value for that grid point. The correction factors are applied to all grid points before the next pass is made.
Usage Notes
1) The oacres function can be quite slow to execute, depending on grid and station data density.
2) The Cressman Analysis scheme can be unstable if the grid density is substantially higher than the station data density (ie, far more grid points than station data points). In such cases, the analysis can produce extrema in the grid values that are not realistic. It is thus suggested that you examine the results of oacres and compare them to the station data to insure they meet your needs.
3) In general, objective analysis is a complex topic, and many schemes for doing it have been developed over the years. The oacres function is provided more as a quick-look feature rather than a rigorous analysis scheme. If you have specific analysis requirements, consider doing your objective analysis outside of GrADS with special purpose programs.
Examples
1) In the simplest case:
oacres(ts,ts.2)
2) To specify your own radii of influence:
oacres(ts,ts.2,12,8,5,4,3,2,1)
pow(expr1,expr2)
The pow function raises the values of expr1 to the power of expr2. No error checking is performed for invalid values in the operands. This function works on both gridded and station data.
Examples
1) To square some value:
pow(expr,2)
2) To duplicate the operation of the mag function:
sqrt(pow(u,2)+pow(v,2))
sin(expr)
Takes the sin of the provided expression. It is assumed the expression is in radiians. Result values are in the range -1 to 1. This function works on both gridded and station data.
skip(expr,skipx,skipy)
Sets alternating values of the expr to the missing data value. This function is used while displaying wind arrows or barbs to thin the number of arrows or barbs.
expr A grid expression.
skipx Skip factor in the X direction.
skipy Number of grid points to skip in the Y direction.
Examples
1) To display every other grid point in both the X and Y direction:
d skip(u,2,2);v
2) To display every grid point in the X direction, but every 5th grid point in the Y direction:
d skip(u,1,5);v
Note that it is not necessary to use the skip function on both the U and V wind components; it is sufficient to populate only one component with missing data values to suppress the plottingof the wind arrow or barb.
smth9(expr)
Performs a 9 point smoothing to the gridded result of the expr.
Usage Notes
1) The result at each grid point is a weighted average of the grid point plus the 8 surrounding points. The center point receives a wieght of 1.0, the points at each side and above and below receive a weight of 0.5, and corner points receive a weight of 0.3.
2) All 9 points are multiplied by their weights and summed, then divided by the total weight to obtain the smoothed value. Any missing data points are not included in the sum; points beyond the grid boundary are considered to be missing. Thus the final result may be the result of an averaging with less than 9 points.
3) If the gridded data is 1-Dimensional, the result is a 3 point smoothing.
sqrt(expr)
Takes the square root of the result of the expr. This function works on both gridded and station data. Values in expr that are less than zero are set to missing in the result.
stnave(expr,dexpr1,dexpr2<,-m cnt>)
Takes an average of station data over time:
expr A valid GrADS expression that gives a station data result.
dexpr1 A dimension expression giving the starting time for the average.
dexpr2 A dimension expression giving the ending time for the average.
-m cnt Optional minimal data count for the average to be taken. If, in the time series, there are fewer available data points for a particular station than the cnt value, then the result for that station is the missing data value. The default cnt value is 1 (ie, even 1 valid station in a time series of even thousands of points would give a valid result for that station).
Usage Notes
1) The times are looped through based on the time interval of the default file. It is thus very important to set the default file to that of the station data file, or a file with the same time interval, or not all station reports will be included in the average.
2) If there is more than one report per station for a particular time, those reports are averaged equally to arrive at a single value for that time. The final average consists of each report for each time being averaged, with missing times not included in the average.
3) Reports from different times are considered to be for the same station when the station id, the latitude, and the longitude all match exactly.
Examples
1) A typical usage of the stnave function would be:
stnave(ts,t=1,t=20,-m 10)
Here an average is taken over 20 times, and if there are fewer than 10 reports for a station then that station will be missing in the final result.
stnmin(expr,dexpr1,dexpr2<,-m cnt)
Examines a time series of station data and returns the minimum value encountered for each station. Operands and usage are the same as the stnave function; see above.
stnmax(expr,dexpr1,dexpr2<,-m cnt)
Examines a time series of station data and returns the maximum value encountered for each station. Operands and usage are the same as the stnave function; see above.
tan(expr)
Applies the trigonometric tangent function to the expr which is assumed to be in radians. Operates on both gridded and station data.
tloop(expr)
When time is a varying dimension in the dimension environment, the tloop function evaluates the expr at fixed times, then reconstructs the time series to obtain a final result that is time varying. The tloop function is required due to the implementation of the GrADS expression evaluation rules, and the implementation of certain other functions. The tloop function can also improve performance for certain calculations.
The tloop function is provided as a way to obtain time series from functions that themselves are not implemented to be able to operate when time is a varying dimension. See the examples below.
Usage Notes
1) The tloop function loops through time based on the time increment of the default file; it is thus important to have the default file set appropriately.
2) The tloop function and the define command work very similarly. In many cases, the define
command can be used to obtain the same result as using tloop. In fact, the define command can be even more useful along those lines, since it also loops through the Z dimension, in effect creating a zloop function. See the define command for more information.
Examples
1) A typical application of the tloop function is to obtain a time series of areal averages using the aave function. Since the aave function will not work when time is a varying dimension, the useof tloop is required:
set x 1
set y 1
set t 1 31
d tloop(aave(ts,lon=0,lon=360,lat=-90,lat=90))
Note that the dimension environment is set up to reflect the kind of plot desired, namely a line plot where time is the varying dimension. Thus it is necessary to fix the X and Y dimensions; the values of those dimensions in this case are not relevent.
2) The tloop function can be used to smooth in time:
set lon -180 0
set lat 40
set lev 500
set t 3 28
d tloop(ave(z,t-2,t+2))
In this example, we are plotting a time-longitude cross section, where each time is a 5 time period mean centered at that time.
3) If we wanted to display a time-longitude cross section (X and T varying), with the data being averaged over latitude, the 'standard' way to do this might be:
set lon -180 0
set lat 40
set lev 500
set t 1 31
d ave(z,lat=20,lat=40)
This calculation could be fairly time consuming, since to perform the average, a longitude-time section is obtained at each latitude. If the time period is long, then this would be a very inneficient operation, due to the ordering of data in a typical GrADS data set. The tloop function might substantially improve the performance of this calculation:
d tloop(ave(z,lat=20,lat=40))
since the average is then done at each fixed time, and is thus just an average of X varying data over Y. Thus the tloop function here is simply being used to force a different ordering to the calculation, although the result is the same.
tvrh2q(tvexpr,rhexpr)
Given virtual temperature and relative humidty, tvrh2q returns specific humidity, q, in g/g. Specifically:
tvexpr A valid GrADS expression where the result is virtual temperature, in Kelvin.
rhexpr A GrADS expression that results in relative humidty, in percent (values from 0 to 100).
This function works only on gridded data.
Usage Notes
1) The conversion formula requires a pressure in mb. tvrh2q assumes that the Z coordinate system is pressure in mb. If Z is a varying dimension, the pressure valid at each grid point is used. When Z is a fixed dimension, the Z value from the current dimension environment is used.
Note that it is possible to provide values from an incorrect pressure level by overriding the current dimension environment:
set lev 500
d tvrh2q(tv(lev=850),rh(lev=850))
In this case, the tvrh2q function would assume a pressure of 500mb, which is the current dimension environment setting for the Z dimension. However, we are providing data from the 850mb level, so the function will produce incorrect results.
tvrh2t(tvexpr,rhexpr)
Given virtual temperature and relative humidity, tvrh2t returns the temperature in degrees Kelvin.
The operation of this function is the same as tvrh2q; refer to the above description for more information.
vint(expr,psexpr,top)
Performs a mass-weighted vertical integral in mb pressure coordinates, where:
expr A GrADS expression for the quantity to be integrated.
psexpr An expression yeilding the surface pressure, in mb, which will be used to bound the integration on the bottom.
top A constant, giving the bounding top pressure, in mb. This value cannot be provided as an expression.
The calculation is a sum of the mass-weighted layers:
where the bounds are the surface pressure on the bottom and the indicated top value on the top. The summation is done for each layer present that is between the bounds. The layers are determined by the different levels of the Z dimension from the default file. Each layer is considered to be from the midpoints between the levels actually present, and is assumed to have the same value throughout the layer, namely the value of the gridpoint at the middle of the layer.
Usage Notes
1) Since the summation is done using the Z levels from the default file, it is important that the default file have the same Z dimension coordinates as the expr.
2) Levels of data below and above the bounds of the summation (surface pressure on bottom; top value on the top) are ignored, even if present.
3) It is assumed the world dimension values for the Z dimension are mb pressure values. The units of g are such that when the expression integrated is specific humidity (q) in units of g/g, the result is g of water per square meter, or essentially precipitable water in mm.
4) It is usually a good idea to make the top pressure value to be at the top of a layer. For example, if the default file (and the data) have pressure levels of ...,500,400,300,250,... then a good value for top might be 275, the value at the top of the layer that extends from 350 to 275 mb.
5) The vint function operates only in an X-Y varying dimension environment.
Examples
1) A typical use of vint might be:
vint(q,ps,275)
to integrate specific humidity to obtain precipitable water, in mm.
If you would like to help test this feature by implementing your own function, please contact doty@cola.umd.edu. Some possible user defined functions might be:
- filtering functions
- grid interpolation functions
- thermodynamic functions